#include "MathematicalUtilities.h"

namespace cagd
{
    double Simpson1D(RealFunction f, double u_min, double u_max, unsigned int div_point_count)
    {
        if (!f)
        {
            return 0.0;
        }

        double step = (u_max - u_min) / (2 * div_point_count);

        double integral = 0.0;

        double even = 0.0;
        for (unsigned int j = 1; j <= div_point_count - 1; ++j)
        {
            double u = u_min + 2 * j * step;
            even += f(u);
        }
        even *= 2.0;
        integral += even;

        double odd = 0.0;
        for (unsigned int j = 1; j <= div_point_count; ++j)
        {
            double u = u_min + (2 * j - 1) * step;

            odd += f(u);
        }
        odd *= 4.0;

        integral += odd;

        integral += f(u_min);

        integral += f(u_max);

        integral *= step;
        integral /= 3.0;

        return integral;
    }

    DCoordinate3 Simpson1D(Function f, DCoordinate3 u_min, DCoordinate3 u_max, unsigned int div_point_count)
    {
        if (!f)
        {
            return DCoordinate3(0,0,0);
        }

        DCoordinate3 step = (u_max - u_min) / (2 * div_point_count);

        DCoordinate3 integral(0,0,0);

        DCoordinate3 even(0,0,0);
        for (unsigned int j = 1; j <= div_point_count - 1; ++j)
        {
            DCoordinate3 u = u_min + 2 * j * step;
            even += f(u);
        }
        even *= 2.0;
        integral += even;

        DCoordinate3 odd(0,0,0);
        for (unsigned int j = 1; j <= div_point_count; ++j)
        {
            DCoordinate3 u = u_min + (2 * j - 1) * step;

            odd += f(u);
        }
        odd *= 4.0;

        integral += odd;

        integral += f(u_min);

        integral += f(u_max);

        integral ^= step;
        integral /= 3.0;

        return integral;
    }
}
